% Sensitivity analysis of delta (discount factor)

  clear;clc;
  close all

% Enter model parameters
  alpha =  3.0;                  % growth function parameter
  beta  =  1.0;                  % growth function parameter
  gamma =  0.5;                  % demand function parameter
  cost  =  0.2;                  % marginal cost of harvest

% Construct state space
  n    = 200;                    % number of states
  smin = 0;                      % minimum state
  sbar = 8;                      % maximum state
  S    = nodeunif(n,smin,sbar);  % vector of states

% Construct action space
  m    = 100;                    % number of actions
  xmin = 0;                      % minimum action
  xmax = 6;                      % maximum action
  X    = nodeunif(m,xmin,xmax);  % vector of actions
  
% Sensitivity analysis of delta
  delta =  0.1:0.1:0.9;                  % discount factor
  
% Empty matrices to store values from change values of delta
  vs    = zeros(n, numel(delta));
  xs     = zeros(n, numel(delta));
  pstars = zeros(n, n, numel(delta));

% Construct reward function
  f = zeros(n,m);
  for k=1:m
     f(:,k) = (X(k).^(1-gamma))/(1-gamma)-cost*X(k);
     f(S<X(k),k) = -inf;
  end

% Construct state transition function
  g = zeros(n,m);
  for i=1:n
  for k=1:m
    snext  = alpha*(S(i)-X(k)) - 0.5*beta*(S(i)-X(k)).^2;
    g(i,k) = getindex(snext,S);
  end
  end

% Pack model structure
  for i=1:numel(delta)
  clear model
  model.reward     = f;
  model.transfunc  = g;
  model.discount   = delta(i);

% Solve infinite-horizon model using policy iteration
  [v,x,pstar] = ddpsolve(model);
  vs (:, i) = v;
  xs (:, i) = X(x);
  pstars (:,:,i) = pstar;
  end

%

plot(S,xs);
title('Optimal Harvest Policy');
legend('Location',"eastoutside");
xlabel('Stock'); ylabel('Harvest');
legend("Box","off");
for i = 1:numel(delta)
    legendInfo{i} = ['\delta = ' num2str(i/10)]; 
end
legend(legendInfo)

%

plot(S,vs); 
title('Optimal Value Function');
xlabel('Stock'); ylabel('Value');
lgd1 = legend('Location',"eastoutside");
legend("Box","off");
legend(legendInfo)

%

sinit =  max(S); 
nyrs = 20;
spaths = zeros(nyrs+1, numel(delta));
for z=1:numel(delta)
spath(:,:,z) = ddpsimul(pstars(:,:,z),n,nyrs);
spaths(:,z) = S(spath(:,:,z));
end

plot(0:nyrs,spaths)
title('Optimal State Path')
xlabel('Year'); ylabel('Stock');
lgd1 = legend('Location',"eastoutside");
legend("Box","off");
legend(legendInfo)

%
% Sensitivity analysis of alpha (growth factor)
%

  clear;clc;
  close all

% Enter model parameters               
  beta  =  1.0;                  % growth function parameter
  delta = 0.9;                   % discount factor
  gamma =  0.5;                  % demand function parameter
  cost  =  0.2;                  % marginal cost of harvest

% Construct state space
  n    = 200;                    % number of states
  smin = 0;                      % minimum state
  sbar = 8;                      % maximum state
  S    = nodeunif(n,smin,sbar);  % vector of states

% Construct action space
  m    = 100;                    % number of actions
  xmin = 0;                      % minimum action
  xmax = 6;                      % maximum action
  X    = nodeunif(m,xmin,xmax);  % vector of actions
  
% Sensitivity analysis of delta
  alpha =  3:0.2:4;                  % growth function parameter
  
% Empty matrices to store values from change values of alpha
  vs    = zeros(n, numel(alpha));
  xs     = zeros(n, numel(alpha));
  pstars = zeros(n, n, numel(alpha));

% Construct reward function
  f = zeros(n,m);
  for k=1:m
     f(:,k) = (X(k).^(1-gamma))/(1-gamma)-cost*X(k);
     f(S<X(k),k) = -inf;
  end

  % Construct state transition function

% Pack model structure
  for l=1:numel(alpha)
  g = zeros(n,m);
  for i=1:n
  for k=1:m
    snext  = alpha(l)*(S(i)-X(k)) - 0.5*beta*(S(i)-X(k)).^2;
    g(i,k) = getindex(snext,S);
  end
  end
  clear model
  model.reward     = f;
  model.transfunc  = g;
  model.discount   = delta;

% Solve infinite-horizon model using policy iteration
  [v,x,pstar] = ddpsolve(model);
  vs (:, l) = v;
  xs (:, l) = X(x);
  pstars (:,:,l) = pstar;
  end

%

plot(S,xs);
title('Optimal Harvest Policy');
legend('Location',"eastoutside");
xlabel('Stock'); ylabel('Harvest');
legend("Box","off");
for i = 1:numel(alpha)
    legendInfo{i} = ['\alpha = ' num2str(3+(i-1)*0.2)]; 
end
legend(legendInfo)

%

plot(S,vs); 
title('Optimal Value Function');
xlabel('Stock'); ylabel('Value');
lgd1 = legend('Location',"eastoutside");
legend("Box","off");
legend(legendInfo)

%

sinit =  max(S); 
nyrs = 20;
spaths = zeros(nyrs+1, numel(alpha));
for z=1:numel(alpha)
spath(:,:,z) = ddpsimul(pstars(:,:,z),n,nyrs);
spaths(:,z) = S(spath(:,:,z));
end

plot(0:nyrs,spaths)
title('Optimal State Path')
xlabel('Year'); ylabel('Stock');
lgd1 = legend('Location',"eastoutside");
legend("Box","off");
legend(legendInfo)